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ABSTRACT 

Aims. We study the dynamics of gas and dust in a protoplanetary disk in the presence of embedded planets. We investigate the conditions for 
dust-gap formation in terms of particle size and planetary mass. We also monitor the amount of dust that is accreted by the planet relative to 
the amount of gas, which is an important parameter in determining the enrichment of solids in giant planets compared to the solid content of 
the central star. 

Methods. We use a new two-fluid hydrodynamics code to solve the flow equations for both gas and dust. For the gas, we use a Godunov-type 
scheme with an approximate Riemann solver (the Roe solver). The dust is treated as a pressureless fluid by essentially the same numerical 
method as is used for the gas. 

Results. We find that it only takes a planet of 0.05 Jupiter masses to open up a gap in a disk with a significant population of mm-sized 
particles. Dust particles larger than 150 fim participate in gap formation. We also find that the formation of the gap severely slows down dust 
accretion compared to that in the gas. Therefore, it is not possible to enrich a newly formed giant planet in solids, if these solids are contained 
in particles with sizes from 150 /im to approximately 10 cm. 



i-H Key words, hydrodynamics - methods:numerical - stars:planetary systems 
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1. Introduction 

Dust plays a major role in planet formation, theoretically as 
well as observationally. First of all, terrestial planets consist 
mainly of solid material. They are giant agglomerates of inter- 
stellar dust particles, and therefore the behaviour of these small 
particles is important during the first stages of planet forma- 
tion. Whether they grow to kilometer-sized bodies by slow co- 
agulation or through the gravitati onal instability of a thin d ust 
layer in the midplane of the disk ( iGoldreich & War3l 19731) is 
still an open question. A mixture of equa l amounts of gas and 
dust i s subject to streaming instabilities ( Yo udin & Goodrnai] 



120051) . and it is not clear whether the resulting turbulent 
struct ure will prohibit grav itational instabilities or enhance 
them jjohansen et al.l l2006'). To answer this question, multi- 
component simulations are needed in which gas and dust are 
able to evolve separately. 

Secondly, giant pla nets can be fo rmed either by gas accre- 
tion onto a solid core jPoUack et alj 1996) or by direct gravi- 
tational fragmentation of the disk ( IBossII1997.) . In the core ac- 
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cretion model, dust growth is of critical importance to build up 
a massive core before the gas disk dissipates. The gravitational 
instability model does not rely on the existence of such a core, 
but still, accretion of dust in a later stage will determine the en- 
richment in solids with respect to t he solar nebula, wh ich is an 
important observational constraint jGuillot et alJ2004i) . 

Finally, dust is the major source of continuum radiation 
in the disk. Therefore, infra-red and sub-millimeter emissions 
mainly indicate the density and temperature of the dust disk. 
The relation of the dust to the density of the gas disk can be 
quite complicated if the disk has a non-trivial structure, espe- 
cially if the dust particles are only marginally coupled to the 
gas. 

It is therefore very important to study dust dynamics 
in gas disks. Axisymmetric disks can be handled analyti- 
cally, and v arious authors have studied dust migration i n 
these disks ('Weid enschiUind Il977t iTakeuchi & Lit] l2002t) . 
Particles move both inward and outward under the influence 
of gas drag, depending on their distai ice from the midplane , 
but the average migration is inward llTakeuchi&LiiJl2002h . 
iHaghighipour & BossI l l2003h studied dust migration in a non- 
uniform axisymmetric nebula, and iTakeuchi & Artvmowic j 
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( 1200 11) investigated the effect of a sharp disk edge on dust mi- 
gration. The overall conclusion is that the structure of the gas 
disk can have a significant effect on the dust distribution, and 
that the distribution of gas and dust can differ substantially. 

When the gas distribution become s more complicat e d, on e 
has to rely on numerical simulations. ISuttner & Yo rke (2001) 
studied dust evolution during the initial collapse of the nebula 
with a two-dimensional multi-fluid h ydrodynamics code to ob- 
tain the large-scale dust distribution. Ijohansen & Kl ahr (2005) 
studied the effect of MHD turbulence on the dust distribution. 
These multi-fluid calculations are computationally very expen- 
sive, requiring more than twice the computing time for ordinary 
gas-dynamical simulations, even if one solves for only a single 
dust particle size. 

To compare dust observations with gas-dynamical simula- 
tions, but at the same time avoid the large computational effort 
of these multi-fluid calculations, one can assume a fixed dust- 
to-gas ratio and solve only for the gas dynamics only. The ef- 
fects of an embedded planet on a gaseous disk have been stud- 
ied numerica lly before by multiple authors (e.g. .Kiev 1999; 
Nelson et^llb oOO: I.u bow et alJ ll99i lO' Angelo et alJ l2002). 
Wolf et a l. (^2002) used these results to simulate observations of 



gaps in protoplanetary disks. Their calculations are valid only 
for the smallest particles that couple extremely well to the gas. 

However, in a protoplanetary di sk where grain growth 
is ver y important (see, for example, IPuUemond & DominikI 
l2005l) . a significant fraction of the dust population is 
large enough to be only mar ginally coupled to the gas. 
IPaardekooper & Mellemal(l2004l) showed that including the full 
gas and dust dynamics leads to a dramatic evolution of the dust 
component in the disk. Planets not massive enough to open 
up a gap in the gas disk do so in the dust disk, which should 
make the orbit of even 0.05 Mj planets visible to ALMA (the 
Atacama Large Millimeter Array). 

In this paper we continue our investigation of dust flow in 
gas disks around embedded planets, wit h mor e attention to de- 
tail flian in [Paardekooper & Mellem^ ll2004 . In Sect. |2l we 
review the equations governing dust flow and gas-dust inter- 
action. Section |3l is devoted to the numerical method, and we 
discuss the initial and boundary conditions in Sect. 0] In Sect. 
|5l we present the results. Section|6lis reserved for a short dis- 
cussion, and we conclude in Sect.0 



2. Basic equations 

Protoplanetary disks are fairly thin, i.e. the vertical thickness 
H is small compared with the distance r from the centre of 
the disk. Typically we use H/r = h = 0.05. It is therefore 
convenient to average the equations of motion vertically and to 
work with vertically averaged quantities only. 

The governing equations are then solved in a cylindrical co- 
ordinate frame (r, (/>), centred on the central star and corotating 
with the embedded planet that has a Keplerian angular velocity 

n. 



2.1. Gas 

The equ ations that control the gas e voluti on are described in 
detail in IPaardekooper & Mellemal (l2005l) . We do not solve 
the energy equation, but instead use an isothermal equation of 
state: 



P = cf E 



(1) 



where Eg is the surface density of the gas and P is the vertically 
averaged pressure. 

The isothermal sound speed Cg is, in hydrostatic equiUb- 
rium, directly related to the disk thickness: 



h VK, 



(2) 



where wk is the Keplerian velocity. 

We use a constant kinematic viscosity which is set by 
assuming an a-parameter (Shakura &Sunvaevll973h of 0.004 
at the location of the planet: 



2.2. Dust 



(3) 



We treat the dust as a pressureless fluid; its evolution is gov- 
erned by conservation of mass and, in absence of external 
sources, conservation of radial and angular momentum. These 
can be written in the following compact form: 



dW dF 

dt dr 



dG 



(4) 



where W is called the state vector and where F and G are the 
fluxes in the radial and azimuthal direction, respectively. S is 
called the source term. 

The components of the state, fluxes, and source-term vec- 
tors can be written in the following form: 



W = r(E. Eu,., Ew^) 
F — r(Eu,., Et;^, Ewr^^) 
G = r(Eu0, EvrV^, Eu^) 



Er2(t;^ + r!)2-Er^ 
-2i:vr{n + V4,) - I 



(5) 
(6) 
(7) 



(8) 



Here, E denotes the dust surface density, Vr and the radial 
and angular velocity, respectively, fi the angular velocity of the 
coordinate frame (corotating with the planet), and $ the grav- 
itational potential of the central star and the planet. The gravi- 
tational potential of the planet is softened over a fraction 0.2 of 
the Roche lobe of the planet. The potential also contains terms 
due to the non-inertial nature of the coordinate system. The 
drag forces are incorporated through f^.r and /d,^. 

The major difference between the equations for the gas and 
the dust is the absence of pressure in the latter In this sense, 
the dust fluid behaves like a gas that is always moving with 
supersonic velocity. This implies that near shock waves, where 
the gas goes from sonic to supersonic flow and where large 
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Fig. 1. Solution to the Riemann problem with initial conditions 
PL — 1-0, ph — 0.5, ul = 0.5, and un — 0.25. The solid 
line represents the dust density and the dots represent the gas 
density, both at t = 0.5. For the gas, the speed of sound was set 
to 1.0. 



density and velocity gradients are present, dust and gas behave 
in very different ways. 

The basic building block for Godunov-type hydrodynamic 
solvers is the Riemann problem, consisting of two stationary 
states separated by a discontinuity. If we consider the two states 
as belonging to neighbouring grid cells, with the discontinuity 
at the interface between the cells, the analytic solution of the 
Riemann problem can be used to define an interface flux be- 
tween the cells. Popular schemes of this ty pe are the Piecewise 
Parabolic Method (PP M, .Woodward & Colella..l984> and the 
Roe solver llRoell98lh . 

For pressureless fluids, the solution to the Riemann prob- 
lem differs significantly from the case with pressure. Figure 
n] shows the analytic solutions for the gas (dots) as well as 
the dust density (solid line), with initial conditions p^ = 1.0, 
PR = 0.5, ul ~ 0.5, and ur = 0.25. The gas solution consists 
of two waves moving at speeds u — and ii + Cg, where 



PL UL 



/PR Ur 



(9) 



If we now let Cg to obtain a pressureless fluid, the sep- 
aration between the waves becomes smaller and smaller, until 
they overlap at Cg = 0. The overlapping waves create a delta- 
function singularity that moves with velocity u (e.g. Leveau3 
[2004"). The singularity arises because Eq.^is only weakly hy- 
perbolic, in contrast with the gas hydrodynamic al equations, 
which are strictly hyperbolic. 

Figure [0 deals with the case ul > ur. When ul < ur, 
the solution consists of the two states moving away from each 
other with a vacuum in between. Note that for both cases, the 



flux is well defined everywhere but at the position of the delta 
function. 

The solution to the Riemann problem not only forms the 
basis for our numerical method (see Sec.|3}, but from Fig.^it 
is also clear that gas and dust react very differently to physi- 
cal shocks, which may lead to gas-dust separation when shock 
waves are present in the computational domain. 

2.3. Gas-dust interaction 

The interaction between the gas fluid and the dust fluid occurs 
only through the drag forces. The nature of the drag force de- 
pends on the size of the particles with respect to the mean free 
path of the gas molecules. We consider only spherical parti- 
cles with radius s. Because the main constituent of the gas is 
molecular hydrogen, for the mean free path we can write: 



A = 



(10) 



When the particle size is small, compared to A, we are in the 
Epstein regime, whereas s ^ A corresponds to the Stokes 
regime. For particles that move subsonically through the gas, 
the transition occurs around A/s = 4/9. For typical midplane 
densities, particles larger than approximately 50 cm move in 
the Stokes regime. In this paper we will focus on particles 
smaller than ^ 1 cm, and therefore we can safely use the 
Epstein drag law. 

In the Epstein regime, when the relative velocity of gas and 
dust is much smaller than the local sound speed, the drag force 
is written as: 



/d 



Ts ' 



(11) 



where 57k is flie Keplerian angular velocity, Av is the velocity 
difference between gas and dust, and Tg is the dimensionless 
stopping time ( Takeuchi & Lin.2002i) : 



^ ^ nppSVK 



(12) 



Here, pp is the particle internal density, s is the radius of the 
particle, and Cg is the sound speed. We have used pp — 1.25 
g cm^'^, and the isothermal sound speed is given by Eq.|2] The 
gas density pg is found by using pg — Sg/2i?. 

For particles that move very supersonically through the gas, 
the drag force becomes: 



•'^ 4V 8 Tg cg 



(13) 



A standar d way to co nnect the supersonic to the subsonic 
regime is (lKwokll975h : 



(14) 



See lWoitke & Heflind ( l2003l) for the general case. 

Only very massive planets of approximately 1 Mj are able 
to accelerate the larger particles to velocities that are compara- 
ble to the speed of sound. For the low-mass planets considered 
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in lPaardekooper & Mellemal ( l2004l) . the subsonic version of the 
drag force proved to be sufficient. 

Particles with Tg ^ 1 are well-coupled to the gas, and the 
radial drift velocity can be written as (.Takeuchi & LirL20 02): 



VrArilt = -VTsVK, 



(15) 



where r/ is the ratio of the gas pressure gradient to the stellar 
gravity in the radial direction: 



i91ogr 



1 



91ogr 



(16) 



The direction in which the dust moves with respect to the gas 
is determined by the sign of rj, and in general the pressure gra- 
dient is negative, making 77 positive, which leads to an inward 
migration of the dust particles. However, if there is a large pos- 
itive density gradient (9 log Eg/9 log r > 1, see Eg. I16>. as 
there is near the outer edge of a density gap, the dust particles 
can move outward. On the other hand, if there is an extra large 
negative density gradient, as is the case near the inner edge of 
a density gap, the dust moves inward with higher velocity. 

Thus, we expect that near the gap edges, the larger particles 
decouple from the gas, altering the dust-to-gas ratio in that re- 
gion. Furthermore, the drag forces may not be large enough to 
prevent gas-dust separation by shocks. Both processes should 
play a role in the dust dynamics near the embedded planet. 

3. Numerical method 

3. 1. Gas and dust advection 

For the evolution of the gas component, we use the R ODEO 
method described bv Paardekooper & Mellenial (l2005l) . It is a 
second-order Eulerian hydrodynamics code, which uses an ap- 
proximate Riemann solver (lRoelll98ll). It derives fro m a gen- 
eral relativistic method ( lEulderink & Mellemal 1 9951) . and it is 
especially suitable to treat non-Cartesian, non-inertial coordi- 
nate systems. A module for Adaptive Mesh Refinement (AMR) 
can be used to obtain high resolution near the planet. 

For the dust, we use the method outlined in iLeveauel(l2004 . 
which has the advantage of being similar to the method used 
for the gas. Below, we only consider the radial direction. The 
azimuthal direction is handled in the same way. 

For every cell interface, we consider the Riemann problem 
defined by a left state Wi^ and a right state Wn with corre- 
sponding fluxes. When Wr.L < < Vr.R, the solution to the 
Riemann problem gives us a vacuum state at the cell interface, 
so we set the interface flux i^int = 0. Otherwise, we compute 
the shock speed v from Eq. |9l and compute the interface flux 
according to 



(Fl + Fr)/2 



if u > 0; 
if u = 0; 
if u< 0. 



(17) 



In the first and the last case, the cell interface is respectively 
left and right of the delta shock, allowing for an easy definition 
of the interface flux. When v — 0, however, the delta shock is 
exactly at the cell interface. We then choose to distribute this 
shock equally over both cells. 



To put this in the language of IPaardekooper & Mellemal 
(I20n5l) . we write F = AW, where the matrix A is given by: 



Vr 
A^ \ Vr 
Vr 



(18) 



The eigenvalues of this matrix are: 

Al = A2 = A3 = Vr. (19) 

The corresponding eigenvectors read: 

ei = (1,0,0) 

62 = (0,1,0) (20) 

63 = (0,0,1). 

A vector A = (Ap, A^, A^) can be projected on these eigen- 
vectors using the following projection coefficients, found by 
solving the system Cb = A, where C is the matrix with the 
eigenvectors: 



h = Ap 
b2 = A^ 

&3 - Arf,. 



(21) 



To approximate the matrix A at an interface of two grid cells, 
we use Roe-averages. If we then project the flux difference 
across the interface we are considering onto the eigenvectors 
of A: 



(22) 



we can then define the first order interface flux just as it was 
defined in the gas case: 



1 



2 ^ 



CTk&kGk, 



(23) 



where dk = sign(Ak). When Ak = 0, we take dk ~ 0. It is 
easy to show that this version of the interface flux corresponds 
exactly with Eq.[nj Having the interface flux in this form, we 
can immediately write down a second-order interface flux that 
is again analogous to the gas case: 



int 



- Fr) - i ^(o-kflk - (c^k - i'k)V'k)Akek, (24) 



where i^k ~ Ak At/ Ar and iJj is the flux limiter . We h ave used 
the same limiter as in Paardek ooper & Mellemal ( l2005l) for both 
gas and dust in all simulations. 

3.2. Accretion 

We model accretion onto the planet by taking mass away from 
the grid at the location of the planet. There are two parameters 
governing the amount of mass taken away at each time-step: 
the accretion area and the accretion rate. We take the accretion 
area to be a circle with a radius equal to a tenth of the Roche 
lobe of the planet, and the accretion rate to be such that the 
time scale for emptying this region is equal to three fifths of 
the orbi tal time scale of t he pla net. The same parameters were 
used in iD'Angelo et alJ (l2002l) and IPaardekooper & MeUemal 

aooi. 
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Fig. 2. Numerical results for the Riemann problem of Fig.[T]for gas (dots) and dust (solid line). The top row shows the density, 
the bottom row the velocity. Left panels: no gas-dust interaction. Middle panels: Tg = 0.5. Right panels: = 0.05. 



3.3. Source terms 

All the source terms exc ept the drag forces are integrated using 
stationary extrapolation llPaardekooper & Mellemal2 005'). The 
combination of exact extrapolation for the angular momentum 
and approximate extrapolation for the other geometrical source 
terms was shown to give reliable results for single-fluid calcu- 
lations of the planet-disk problem. Because these results were 
not affected by increasing the resolution, we can conclude that 
using this treatment of the source terms is accurate enough. 
However, when other forces besides gravity are included, it 
is not guaranteed that a stationary solution exists. Simis et al. 
( I2OOII) found that the source terms arising in two-fluid calcula- 
tions are best integrated analytically. Therefore, the drag forces 
are incorporated separately using an ordinary differential equa- 
tion: 



(25) 



During this step, we keep the dust and gas density fixed. Using 
Eq. ^5 we can write an ordinary differential equation for the 
velocity difference between gas and dust: 



dt 



(Av) 



-a 



v/l + 13 Av^Av, 



(26) 



where a = ^ and /? = j 
analytically to give 



28c2 



. This equation can be integrated 



Av{t) 



2Ava exp {—at) 



(27) 



Note that when /3 = 0, we recover the ordinary exponential 
decay of the velocity difference. Using conservation of total 
momentum: 



Si; 



^gVg 



constant, 



(28) 



we can solve for the evolution of the velocities of gas and dust 
separately, while keeping the total momentum conserved ex- 
actly. 

3.4. Test problem 

Gasdynamical codes are often tested against a simple Riemann 
problem ( Sod 1978) . also known as the shock tube. In this sec- 
tion, we test the method for dust advection against the Riemann 
problem shown in Fig.[n using a 100-zone grid. We also study 
the effect of the interaction between gas and dust. The ini- 
tial dust-to-gas mass ratio is 1 . The simulations were run until 
t — 0.5, which took approximately 50 time-steps. 

In Fig. 12] we show the results for three different stopping 
times: Tg = 00 (left panels), Tg = 0.5 (middle panels), and 
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Tg ~ 0.05 (right panels). We see that for the case of no interac- 
tion, we reproduce the result of Fig. for gas and dust. The 
delta shock is smeared over approximately 4 computational 
cells, but it moves at the correct speed. Note that because the 
velocity is less than the sound speed everywhere, the waves in 
the gas are not shock waves, but are simply sound waves. 

For the middle panels, the stopping time is comparable to 
the sound-crossing time and the dust velocity increases towards 
the gas velocity, while by conservation of momentum, the gas 
velocity decreases with respect to the left panels. The veloci- 
ties approach each other the fastest at the position of the delta 
shock which has decreased in amplitude because the gas and 
dust densities also approach each other 

When we decrease the stopping time by a factor of ten, the 
densities and velocities of gas and dust are almost equal every- 
where. In the right panels of Fig.|2] there is only a small bump 
left of the delta shock, and the gas has lost enough momentum 
to the dust to visibly slow down the sound waves. Note that this 
only happens because the gas and dust densities were equal ini- 
tially. For an interstellar dust-to-gas ratio of 1 : 100, the gas is 
not affected by the dust. 

3.5. Limits 

The approach presented here cannot be used in all circum- 
stances. To begin with, the continuum form for the equations 
of motion is only valid if two requirements are met. First, there 
need to be enough particles in the volume under consideration 
(in our case, one grid cell) to work with average states. Second, 
there need to be enough interparticle collisions to make such 
an average state meaningful. For a dust fluid, both criteria 
are only fullfilled if the particles are well-coupled to the gas 
jGaraud et al. 1'2004'). In other words, we can only work with 
particles for which Tg <C 1 everywhere in the computational 
domain. For larger dust boulders, a particle-based method 
is mo re appropriate ilohansen et alJl2005t iFromang & Nelsorj 
1200,4 . 

On the other hand, when the stopping time is very small, 
the source terms due to gas-dust friction become stiff. Even 
though we integrate these source terms in an implicit way, 
which makes the method stable for these small particles, it does 
not necessarily give us the correct answer. This is because the 
gas source terms due to pressure are dealt with in the advection 
step, and pressure is precisely what causes gas-dust separation. 
Therefore, even for Tg = 0, the dust would numerically slowly 
drift inward, while according to Eq. its radial velocity 
should be zero. The range in Tg that we can consider is there- 
fore also bounded from below. The time scale for the advection 
(the hydrodynamical time-step, determined by the Courant- 
Friedrichs-Lewy (CFL) condition) should not be much larger 
than the time scale for gas-dust coupling. Combining the two 
limits, we can conclude that we can only consider particles for 
which the stopping time satisfies: 

17k Ai « < 1. (29) 

In practice, for the gas densities considered here (see Sect. |^, 
Eq. |29] limits the particle size from below at approximately 



150 /im and from above at approximately 1 cm. In all our sim- 
ulations we checked that lowering the time-step At did not in- 
fluence the results. 

It is important to realize that Eq.|29lshould hold true every- 
where on the grid. When the gas density changes significantly, 
for example due to the formation of a gap, the stopping time can 
vary by orders of magnitude. The condition Tg <C 1 for gap- 
opening planets makes the simulations very expensive. To do 
the calculations in a reasonable amount of time, we restrict our 
parameter space to non-gap-opening planets. Note, however, 
that gap-opening planets also pose problems for particle-based 
methods because very many particles are needed to achieve 
enough resolution inside the gap. 

The two-dimensional nature of these simulations also 
puts conditions on the size of the parameter space we can 
consider When the Roche lobe of the planet is smaller 
than the disk thickness, three-dimensional effects come into 
play, lowering accreti on and migration rates for example, 
(D'Angelo et all l2003h . However, global three-dimensional 
multi-fluid calculations are beyond current computational re- 
sources. Furthermore, because of dust settling, the dust disk 
will be much thinner than the gas disk, depending on the par- 
ticle size. The two-dimensional approach will therefore cause 
fewer discrepancies in the dust distribution. Still, one has to be 
very careful when interpreting results for planets with M < 0.1 
Mj. 

4. Initial and boundary conditions 

The standard numerical resolution we use is {nr,nci,) = 
(128, 384). This way, the cells close to the planet are of equal 
size in the radial and in the azimuthal directions. We do not re- 
solve the Roche lobe using a low resolution like this, but this 
resolution is sufficient to capture the long-term global evolu- 
tion of the disk. For the runs with dust accretion, we used 4 
levels of AMR to achieve a resolution that is 16 times higher 
near the planet. 

As mentioned before, we use an aspect ratio of ft. = 0.05, 
and the initial gas surface density distribution is constant, with 
Sg = 34 g cm^^. This surface density is appropriate for a disk 
with 0.01 M0 within 100 AU at approximately 7 AU, between 
the location of Jupiter and Saturn in the solar system, or at 13 
AU in the Minimum Mass Solar Nebula. Note, however, that 
the fundamental parameter is the stopping time, and our results 
can therefore be scaled to disks of arbitrary mass by adopting a 
different grain size for a given Tg. 

The initial dust-to-gas ratio is 0.01. Our distance unit is the 
radius of the planet's orbit, and in this unit the inner boundary 
is at r = 0.4 and the outer boundary at r = 2.5. 

It is easy to see that for a constant surface density the stop- 
ping time is also constant: 

Tg = 0.0461 f —) , (30) 
Vcm/ 

where the size of the particles s is given in cm. 

The gas rotates with a slightly sub-Keplerian velocity due 
to the radial pressure gradient (note that even when there is 
no density gradient, there is a pressure gradient caused by the 
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Fig. 3. Relative radial velocity of gas and dust (1 mm) in a disk 
without a planet. The dots represent the numerical solution, and 
the solid line represents Eg. 1151 



gradient in the sound speed). The dust rotates exactly with the 
Keplerian velocity, initially. All radial velocities are taken to be 
zero. 

The constant kinematic viscosity is set so that a = 0.004 at 
r = 1 in all models. For the disk described above, the critical 



I lanet mass for gap opening is about 1 Jupiter mass (1 Mj) 
Brvden et al. ^1999). We varied the mass of the planet between 
0.001 Mj and 1.0 Mj. 

The boundary con di tions are the same as in 
|Paardekooper & Mellema' ('2005^, deriving from 'Godon| 
il996) . They are specifically designed to be non-reflecting, 
so that all waves generated in the simulated region leave the 
computational domain without influencing the interior 



5. Results 

5. 1. Axisymmetric disk 

To begin with, we consider a disk without a planet. The ana- 
lytic solution for the drift velocity is given in Eq. ^] Figure 
|3] shows this analytic solution for particles of 1 mm, as well 
as the numerical solution after 10 orbits of the gas at r = 1. 
The analytical solution is indistinguishable from the numerical 
result, showing that the method handles the large source terms 
arising due to the small Tg very well. The angular velocity of 
the dust is always equal to the angular velocity of the gas. The 
Coriolis force starts to play a role only when the relative radial 
velocity becomes very large. 



5.2. Dust response to a spiral wave 

An embedded planet perturbs the gas disk by launching spiral 
density waves. In this section, we investigate the dust response 
to such a single wave. In Fig.|5 we show a close-up on part of 
the spiral pattern generated by the planet. The velocity arrows 
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Fig. 4. Close-up on a spiral density wave, with relative veloc- 
ity arrows superimposed. Outside of the spiral wave, the dust 
drift is directed inward, according to Eq. Close to the spi- 
ral, the particles move towards the centre of the wave. For 5 mm 
dust particles, as shown here, the maximum relative velocity is 
2.6 10^^, in units of the orbital velocity of the planet. 



indicate the relative velocity of gas and dust. In the upper right 
direction, the density and velocity of the gas are close to their 
initial values, and in this region, the dust drifts slowly inward, 
according to Eq.^] In the lower left part of the figure the dust 
particles drift inward again, but at a higher speed, due to the 
lower gas density. At the location of the spkal wave, gas density 
and velocity are such that the dust always moves towards the 
centre of the wave. This is not surprising in view of Eq. ^] 
which states that in an axisymmetric disk, dust particles will 
always move towards the highest pressure. Note that even in 
this perturbed disk, the angular velocities of gas and dust are 
still almost equal: the particles move only radially through the 
gas. Therefore, the dust response to a spiral wave is the same 
as the response to an axisymmetric density wave, but stretched 
in the radial direction. The effect of a density wave, as shown 
in Fig.|4] remains very local and is of minor importance for the 
global disk evolution. 

5.3. Global disk evolution 

IPaardekooper & Mellemal (|2004l) showed that planets of ~ 0.1 
Mj are able to decouple the evolution of the gas disk and the 
dust disk. More specifically, they showed that the dust particles 
tend to move away from the orbit of the planet, creating a deep 
annular gap. Furthermore, the final dust disk appeared to be 
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Fig. 5. Surface density of gas and dust (s = 0.1 cm) after 500 
orbits of a 0.1 Mj planet at cf) — (opposite to the planet). 
The dust density is multiplied by 100. The dots indicate the 
positions of the 2:3, the 3:2, and the 2:1 resonances. 

structured near mean motion resonances. In this section, we 
further explain the mechanism for structuring the disk. 

First of all. Fig. |5] shows a radial density cut through a 
di sk perturbed by a 0.1 Mj planet after 500 orbits. Again, as 
in IPaardekooper & Mellemal (l2004l) . we see a clear gap and a 
density enhancement at the position of the 2:1 mean motion 
resonance. 

Note that since Tg ^ 1, the dust particles are never influ- 
enced by the planet directly at the radii of the indicated res- 
onances, they only react to the gas density and velocity. The 
mechanism for structuring the dust disk is therefore intrinsi- 
cally different from ordinary resonance capture, which is the 
usua l explanation for the existence of the Plutinos in the Kuiper 
he1t(lT.mi&.TewitJ2002h. 

In Fig. |6j we plot the radial acceleration of dust particles, 
excluding the drag force: 



or 



(31) 



Figure|6lshows the azimuthally averaged a, with a close-up on 
the outer region of the disk. The first thing to note is that the 
dust is always accelerated inward, except for the region around 
r = 1.1. Inside the planet's orbit, the acceleration is very nega- 
tive, indicating that the dust will quickly move toward the inner 
boundary. 

The large inward and outward accelerations near r = 0.9 
and r = 1.1, respectively, are due to the inner and outer edges 
of the density dip carved in the gas. The pressure gradient at 
the outer edge is large enough to push the dust outward. The 
densest ring outside the planet's orbit forms where a switches 
sign at r = 1.3. Particles approach this location in the disk 




Fig. 6. Azimuthally averaged radial acceleration a of dust par- 
ticles after 100 orbits of a 0.1 Mj planet. The close-up shows a 
in the outer region of the disk. The 2: 1 mean motion resonance 
is visible as a small bump in a. 

from both smaller and larger radii, which leads to a strong dust 
accumulation. 

The close-up in Fig.|6lindicates the outer 2:1 mean motion 
resonance, which shows up as a bump in the radial acceleration 
at r = 1.6. Note that a is still negative near the resonance, so 
the dust keeps moving inward. However, near the resonance, 
the dust particles are decelerated, and therefore the dust piles 
up at that location. From Fig.|5] we can also see that the density 
peaks just outside the gap are close to, but not exactly located at 
the position of the 3:2 and the 2:3 resonance, showing that the 
mechanism by which they are formed is different. In contrast 
to the region close to the shallow gas gap, there is no apparent 
pressure gradient near the 2:1 resonance to structure the dust 
disk. Nonetheless, the resonance does show up in the acceler- 
ation in Fig.|6l and it is therefore the velocity field rather than 
the density of the gas that makes the dust move into the 2:1 
resonance. 

The adopted initial gas surface density is correct for a dis- 
tance of 13 AU in the Minimum Mass Solar Nebula (MMSN). 
For this location of the planet, the initial dust mass between 
r = 1.2 and r = 1.5 is 5.4 Earth Masses (M®). After 500 
orbits of the planet, the dust mass in this region has grown to 
10.9 Mg, an enhancement of more than a factor of 2 (see Fig. 
|5}. This may have serious implications for building planetes- 
imals (boulders with sizes larger than 10 centimeters) in this 
region. 

We have not included dust diffusion in our models. 
Turbulent transport of small dust grains is not yet fully under- 
stoo d, but it can be described b y diffusion for the smallest parti- 
cles djohansen & Klahr"2005). To see what effect diffusion has 
on dust-gap formation, we gave the dust fluid the same turbu- 
lent diffusion coefficient as the gas, and included the viscous 
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force terms in the dust source update. It turned out that dust 
diffusion affects the dust density distribution only at the 1 % 
level, indicating that this is not an important process. However, 
the largest particles we con sider may be s ubject to vortex trap- 
ping in turbulent eddies jjohansen & Klah32005l) . and for these 
particles, the diffusion approach may not be entirely valid. 

5.4. Dependence on particle size 

We have seen that particles of 1 mm decouple from the gas 
due to the disk perturbation induced by a small planet. As we 
examine smaller and smaller particles, we expect a certain min- 
imum particle size for this decoupling to occur because smaller 
particles couple better to the gas. We can make a simple es- 
timate for this Sniin by looking at the region in the disk near 
r = 1.1, where the dust experiences the largest positive accel- 
eration a,„ax- According to Fig.|6j Cmax = 0.005. The radial 
gas velocity in a standard accretion disk is 

vr,g = (32) 

The gas radial velocity in a disk with a planet is similar. 
Balancing the outward acceleration with the radial drag force 
gives an equilibrium velocity difference between gas and dust 
of 

Ts 

AVr = ftinaxTT"- (33) 

When this velocity difference is larger than the radial gas veloc- 
ity, dust particles can be accelerated outward around r = 1.1, 
creating a low density gap and a high density ring. The mini- 
mum stopping time for this to occur is 



2rar, 



(34) 



For our disk parameters, this implies a minimum particle size of 
approximately 500 /im. Note, however, that although smaller 
particles will not be pushed outward, they will travel inward at 
a different speed near the edge of the gas dip. As in the case 
for the 2: 1 resonance, where in spite of the acceleration being 
negative everywhere a dust ring still arises, these particles will 
eventually create a gap, but on a much longer time scale. 

Figure0shows the radial dust density for four different par- 
ticle sizes after 100 orbits. First of all, it is clear that the larger 
particles react even more dramatically on the planet. The three 
bumps at r — 0.7, r = 1.3, and at corotation all grow very fast 
with increasing particle size. After 400 orbits, the dust-to-gas 
ratio of 5 mm particles at r = 1.3 is enhanced by more than an 
order of magnitude. The width of the dust gap does not depend 
on the size of the particles. 

The outer disk (r > 1.1) reaches a quasi-static state after 
approximately 100 orbits. From then on the only evolution is in 
the further growth of existing features. The inner disk is slowly 
cleared, however, but on a time scale that is set by the unper- 
turbed motion of gas and dust. If we take Eq.|32]for the velocity 
of the gas and Eq.^]for the relative velocity, we can write the 
following for the dust velocity: 




,2" 



h?VK- 



(35) 



Fig. 7. Radial density cut opposite to the planet for four differ- 
ent particle sizes after 100 orbits of a 0.1 Mj planet. The dust 
densities are multiplied by 100. 



This gives a time scale for clearing the region inside r = 1 of 
approximately 50000 orbits. 

Looking at the distribution of 0.1 mm particles in Fig.0 we 
see that even these well-coupled particles start to create a clear- 
ing in the inner disk. However, since this planet is not strong 
enough to keep all the dust outside of r = 1.3, the inner disk is 
still fed from outside. Particles larger than 500 ^va. are able to 
create a fast gap for this planet. This agrees with the estimate of 
Eq.|34] Based on Fig.0 we estimate that particles larger than 
150 /im will open up a dust gap within 500 orbits. 

5.5. Dependence on planetary mass 

Smaller planets do not dramatically perturb the disk, so we ex- 
pect that a minimum planetary mass for dust-gap opening also 
exists. In this section, we consider only 1 mm particles. 

Figure|8]shows the radial dust distributions for four differ- 
ent planetary masses. The smallest planet is not able to produce 
a density dip near its orbit, so the gas density is equal to 1 ev- 
erywhere except in the spiral waves. From Fig.|8] we conclude 
that a planet of 0.05 Mj is the lowest-mass planet that is able to 
create a gap, albeit on a longer time scale than for more massive 
planets. 

The highest mass planet in Fig. [Sjis on the edge of cre- 
ating a gap in the gas disk. As a result of the lower gas den- 
sity, the stopping time for particles near the orbit of the planet 
is relatively large, which allows for the strong evolution near 
r = 1. The time scale for dust-gap formation for this planet 
is about half the time scale for gas-gap formation, which is 
approximately 100 planetary orbits. Also, due to the strong 
spiral waves, the 2:1 mean motion resonance plays a major 
role; the resonance is able to suck material up from both sides, 
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Fig. 8. Radial density cut opposite to the planet for four differ- 
ent planetary masses after 100 orbits and for a particle size of 1 
mm. The dust densities are multiplied by 100, and the planetary 
masses are in units of M j 



leading to a less pronounced peak near r = 1.4 and to an 
empty region around r = 1.8. Figure |9l shows the dust sur- 
face density for this planet. Note that the resonance-induced 
gaps at r — 0.55 and r — 1.8 are as deep as the gap near the 
planet's orbit. Observationally, it may be difficult to distinguish 
a single-planet system like this from one with multiple planets. 

Because of the conditions discussed in Sec. 13.51 a 0.5 Mj 
planet is the maximum planetary mass we can accurately model 
with a two-fluid calculation. However, based on the relation be- 
tween the radial acceleration a and the resulting dust distribu- 
tion (compare Figs.|5]and|6|l, we can still make a few remarks 
concerning higher-mass planets. 

First of all, we expect similar things to happen for a 1 Mj 
planet as for a 0.5 Mj planet, only in a more extreme form. 
Secondly, we noted in Sec. 15. 31 that the dust ring just outside 
the gap forms where a changes sign from positive to negative. 
In Fig.^l we show a together with the gas surface density S 
for a 1 Mj planet. The gas gap clearly stands out, and we see 
that a changes sign at the outer edge of the gap around r ^ 1.4. 
A dust ring at that location would coincide with the outer edge 
of the gas gap. Therefore, a possible dust gap would have the 
same width as the gap in the gas disk. 

Also from Fig.[lOl we see that a is structured up to r « 2.5, 
which is larger than for the 0.1 Mj planet (see Fig.|6}. For a 
planet at Neptune's orbit (semi-major axis 30 AU), the resonant 
structure in the dust disk can therefore extend to at least 75 AU. 

5.6. Flow within tlie Roctie lobe 

The flow of dust deep within the Roche lobe (Rr ) of the planet 
is relatively simple. Due to the strong pressure gradient in the 




Fig. 9. Surface density of 1 mm dust particles after 100 orbits 
of a 0.5 Mj planet. 




Fig. 10. Azimuthal averages of the radial acceleration a (solid 
line) and gas surface density E (dashed line, shown for 0.4 < 
r < 1.8) for a 1 Mj planet after 400 orbits. 

direction of the planet, the relative velocity of the dust particles 
is always directed towards the planet. In other words, if we take 
a cylindrical coordinate frame (s, 6) centred on the planet, then 
Avs — Vs,d — Vs.g is always negative. 

In Fig. EI we show this relative velocity as a function of 
distance to a 0.1 Mj planet for 1 mm particles. The first thing to 
note is that the relative velocity changes sign very close to the 
planet. That is because this was a run with gas and dust accre- 
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Fig. 11. Scatter plot of the relative velocity in the direction of 
a 0.1 Mj planet. The distance to the planet is in units of the 
Roche lobe Rr , and the unit of relative velocity is the Kepler 
velocity at the position of the planet. The size of the dust parti- 
cles in this simulation is 2 mm. 



tion, and therefore most of the^as ins ide 0.1 Rr was removed 
(see IPaardekooper & Mellema ( 2005 ) for the accretion proce- 
dure), leading to a reversed pressure gradient locally. Outside 
the accretion area, dust moves slowly towards the planet with 
respect to the gas. The spread in velocities stems from the fact 
that the gas density distribution around th e planet is n ot sym- 
metric with respect to the planet (D'Angelo et alJ2002i) . 



If we take the average relative velocity in the Roche lobe 
as 4.0 10^^ in units of the Kepler velocity at the location of 
the planet, then we get a time scale to clear the Roche lobe of 
dust of less then 100 orbits for particles of 1 mm. This is an 
important time scale for accretion, as we will see in the next 
section. 



5. 7. Gas and dust accretion 

We now turn to the interesting problem of dust accretion. The 
accretion of dust particles compared to the gas is important for 
the final enrichment of gaseous giant planets. When the dust 
accretion rate is higher than 0.01 (the initial dust-to-gas ratio) 
times the gas accretion, the planet will get enriched in solids, 
and when the dust accretion is lower than 0.01 the gas-dust 
mixture that is accreted onto the planet is relatively gas-rich. In 
the discussion below, we multiply the dust content of the disk 
by 100, so that gas and dust accretion can be directly compared. 

The mass range we can consider is limited because for 
Mp < 0.1 M.T , three-dimensional effects start to play an im- 
portant role llD' Angelo et alJl2003b . while planets larger than 
approximately 0.5 Mj open up a deep gas gap in the disk, 
which we cannot treat in the fluid approach. Therefore, to elu- 
dicate the basic principles that govern dust accretion and to pro- 
vide a starting point for future, more detailed simulations, we 
only consider a 0.1 Mj planet. 

In the left panel of Fig.^J we show the accreted mass of 
gas and dust in units of gas and dust disk masses as a function 
of time. The accretion of gas (solid line) w as extensively dis- 
cussed in Paard ekooper & Mellem a (2005) and is shown here 
only for comparison. After approximately 200 orbits, it reaches 
a constant value of 1.5 lO""* gas disk masses per orbit. 

The accretion rates for dust particles of different sizes show 
a different behaviour. In the first few orbits, relatively more 
dust is accreted than gas, which means that the dust mass that 
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is accreted is larger than 0.01 times the gas mass that is ac- 
creted. This is the stage where the immediate surroundings 
of the Roche lobe are cleared of dust, due to the strong pres- 
sure gradients in the planetary atmosphere. Later, the compet- 
ing process of dust-gap formation sets in, and consequently, 
dust accretion slows down. In the end, for all particle sizes, rel- 
atively more gas is accreted than dust, which means that the 
planet will not be enriched in solids of these particular sizes. 

The right panel of Fig.[21shows the ratio of accretion rates 
for the different particle sizes to the gas accretion rate. We see 
that initially the accretion rate of dust is higher than 0.01 times 
the accretion rate of the gas or all particle sizes, correspond- 
ing to the Roche lobe clearing mentioned before. If we look 
at the accretion rate of 1 mm particles, after the first sharp 
peak, the accretion rate stays almost constant for 30 orbits, af- 
ter which the process of gap formation sets in, and the accretion 
rate slowly declines. After 300 orbits, the accretion rate for the 
dust is less than half the accretion rate of the gas, and it is still 
declining. 

For the larger particles, we see a larger initial accretion rate 
because these more weakly coupled dust particles move up the 
strong pressure gradient near the planet more easily. However, 
for larger particles, gap formation sets in earlier, as well, and 
it will eventually win. A particle size of 2 mm results in the 
largest accreted mass of dust in this first stage of accretion. 
For 4 mm particles, gap formation sets in before the planet can 
accrete a significant amount of surrounding dust. However, in 
the long run, more small than large particles will have been 
accreted, as can be seen from the left panel of Fig. 1121 

The planet clears the outer part of the gap of 4 mm par- 
ticles in approximately 100 orbits (see Fig.0. At this point, 
the accretion rate curve levels off at a value of 0.2 Mgas, but 
around 150 orbits, it starts to decrease again. Accretion around 
this time comes from the dense ring at r = 1, which will slowly 
diffuse into the inner disk. After it is gone, dust accretion will 
be negligible compared to gas accretion. 

6. Discussion 

Dust gaps opened by low-mass planets have some in- 
teresting observational a nd the oretical consequences. In 
IPaardekooper & Mellemal ( l2004l) . it was shown that the 
Atacama Large Millimeter Array (ALMA) will provide the res- 
olution to observe these gaps in protoplanetary disks at the dis- 
tance of the Taurus starforming region. Observationally, the rel- 
atively large particles we consider are only important at mm- 
wavelengths. To create gaps and inner holes in disks that can 
be detected at infra-red wavelengths, more massive planets are 
needed to remove the smallest particles. 

Next, we will explore some more possible consequences for 
planet formation. Before we can make any quantitative state- 
ments about the importance of the limited size range we con- 
sider in this paper, we need to determine what fraction of the 
total dust mass these particles represent. If we look at the size 
distribution from Fig. 7 of Dullemond & Dominik (2005) we 
see a bimodal distribution (dust and planetesimals), with a tran- 
sition around particles of 10 cm. It is interesting to note that 
for typical mid-plane densities, a size of 10 cm gives rise to 



Ts « 1 (see Eq. I30> . which is the upper limit for which the 
fluid approach is valid. Boulders larger than 10 cm are not 
influenced by the gas as much as the smaller particles, and 
therefore will not participate in gap formation due to pressure 
gradients, while particles smaller than 150 /im are too well- 
coupled to the gas to form a gap (see Sec. \5.4l . From Fig. 7 
of Dullemond & Dominik (2005), we estimate that half of the 
total mass of particles smaller than 10 cm is in particles larger 
than 150 /im. Furthermore, about 30 % of the total dust mass 
is in particles with sizes between 150 /im and 10 cm. 

The above estimate implies that the inner disk is denied 
30 % of the amount of solids it would receive from the outer 
disk without the presence of the planet. If (terrestrial) planets 
are still being formed in the inner disk, it may slow down this 
process considerably. 

Dust piles up outside the planet's orbit the (see Fig.|5|i. In 
Sect. 15.31 we mentioned that after 500 orbits of the planet, the 
amount of dust had already doubled. Because the planet acts as 
a barrier for particles moving inward, more and more dust will 
stream into the resonance from the outer disk. It is interesting 
to note that this is in the region where the Kuiper Belt resides 
with respect to Neptune. This large amount of dust may well 
have influenced planetesimal formation in that region. 

If an already-formed giant planet is present in the inner 
disk, this planet is also denied 30 % of the amount of solids 
it would be able to accrete without the presence of the outer 
planet. One could imagine that in the present configuration of 
our Solar System, with three dust-gap-opening planets outside 
5.2 AU, the amount of solid material accreted by Jupiter would 
be dramatically small, even without the mechanism outlined in 
Sect. Ell 

Depending on the dust size distribution in the disk, the low 
relative accretion rates for mm-sized particles found in Sect. 
15.71 will have a significant effect on the final composition of 
the planet. Our results show that the accretion rates for these 
particles are, on the long run, more than an order of magnitude 
lower than the gas accretion rate. This will significantly affect 
the final enrichment in solids for giant planets, even more so if 
a dust-gap-opening planet exists in the outer disk. 

In other words, if, as is standard in planet-formation mod- 
els, one considers gas accretion and the accretion of planetes- 
imals separately, the gas that is accreted is relatively poor in 
solids (by a factor of 2) compared to the solar value. This pro- 
cess wo rks against th e usual enrich ment scenari os for the giant 
planets (I pautier etal...2001a^b:.Owen et al.ll999l) . However, as 
was mentioned before, we cannot make firm statements about 
planets of higher mass. More definitive conclusions about dust 
accretion onto giant planets require further study. 

Dust growth may also play an important role in structuring 
the disk. When the growth time scale Tg is long compared to the 
dust-gap-formation time scale, what effectively happens is that 
particles will be removed from the dust gap as soon as they be- 
come large enough. This means that the gap will also be cleared 
of the smallest particles on a time scale Tg. However, the pro- 
cess of dust growth in disks is not understood well enough to 
make any quantitative statements about it. 

It is important to realize that the dust-gap-formation pro- 
cess basically only needs a shallow gas dip formed by an em- 
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bedded planet. As long as more complicated disk models in- 
cluding magnetic fields and radiative transfer do not change 
this gas structure, the basic picture outlined in this paper will 
not change. Three-dimensional effects may be important, how- 
ever, because gas and dust behave differently in the vertical 
direction: dust particles will settle to the midplane, while the 
gas forms a pressure-supported structure. Angular momentum 
transport, and therefore also gap formation, is reduced in three- 
dimensional models, compared to vertically integrated models, 
but because of settling effects, the dust will mainly interact 
with gas near the midplane where the planetary potential is 
strongest. It is not clear how these competing processes will af- 
fect our results; three-dimensional, multi-fluid simulations are 
needed to answer that question. 

7. Summary and conclusions 

In t his paper, we have elaborated on the findings of 
IPaar dekooper & Mellema (2004), who showed that intermedi- 
ate mass planets are able to open dust gaps in gas disks. We 
have shown how the formation of the gap depends on particle 
size and planetary mass, indicating that particles larger than 
150 /xm and planets larger than 0.05 Mj are needed to open a 
gap in a typical protoplanetary disk. Planets larger than 0.5 Mj 
clear a gas gap around their orbit, increasing the stopping time 
Ts locally by 2 orders of magnitude, leaving the dust inside the 
gap largely decoupled from the gas. If the planet clears a gas 
gap, the dust gap will be of equal width. 

Within the Roche lobe of the planet, the pressure gradient 
is such that solids will quickly move to the centre of the plan- 
etary atmosphere, allowing for a very high dust accretion rate. 
However, this enhancement is more than counterbalanced by 
the local depletion of dust due to the formation of the dust gap. 
In the end, much more gas can be accreted compared to solids. 
This may have serious consequences for the enrichment of gi- 
ant planets. 
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